We are IntechOpen, the world's leading publisher of Open Access books Built by scientists, for scientists

Open access books available 5,300

130,000 155M

International authors and editors

Downloads

Our authors are among the

most cited scientists 154 TOP 1%

Selection of our books indexed in the Book Citation Index in Web of Science™ Core Collection (BKCI)

# Interested in publishing with us? Contact book.department@intechopen.com

Numbers displayed above are based on latest data collected. For more information visit www.intechopen.com

# **Numerical Geodynamic Modeling of Continental Convergent Margins**

Zhonghai Li 1,2,3 , Zhiqin Xu <sup>2</sup> and Taras Gerya 3

<sup>1</sup>*FAST Laboratory, CNRS/University of Paris 6 and 11* <sup>2</sup>*State Key Lab of Continental Tectonics and Dynamics, Institute of Geology, Chinese Academy of Geological Sciences* 3 *Institute of Geophysics, ETH-Zurich* <sup>1</sup>*France* <sup>2</sup>*China* <sup>3</sup>*Switzerland*

# **1. Introduction**

The continental convergence (subduction/collision) normally follows the oceanic subduction under the convergent forces of lateral ridge push and/or oceanic slab pull (Turcotte and Schubert, 2002). During these scenarios, a large amount of positively buoyant materials enter the trench causing slow down of the convergence that, eventually, may stop. However, before collision ceases, convergence between the plates can continue actively for tens of millions of years after ocean closure as it is testified by the 50 *Ma* active collisions in the Western Alps and Himalayas (e.g. Yin, 2006).

A remarkable event during the early continental collision is the formation and exhumation of high-pressure to ultra-high-pressure (HP-UHP) metamorphic rocks, which is one of the most provocative findings in the Earth sciences during the past three decades. Occurrences of UHP terranes around the world have been increasingly recognized with more than 20 UHP terranes documented (e.g. Liou et al., 2004), which have repeatedly invigorated the concepts of deep subduction (>100 *km*) and subsequent exhumation of crustal materials (e.g. Chopin, 2003). It has been suggested that the HP-UHP metamorphism can be considered as a "hallmark" for the modern plate tectonics regime characterized by colder subduction and started from a Neoproterozoic time (e.g. Brown, 2006, 2007).

The understanding of the dynamics of continental convergent margins implies several different but strictly correlated processes, such as continental deep subduction, HP-UHP metamorphism, exhumation, continental collision and mountain building. Besides the systematic geological/geophysical studies of the continental convergent zones, numerical modeling becomes a key and efficient tool (e.g. Burov et al., 2001; Yamato et al., 2007; Gerya et al., 2008; Warren et al., 2008a,b; Li and Gerya, 2009; Beaumont et al., 2009; Li et al., 2011). The tectonic styles of continental subduction can be either one-sided (overriding plate does not subduct) or two-sided (both plates subduct together) (Tao and O'Connell, 1992; Pope and Willett, 1998; Faccenda et al., 2008; Warren et al., 2008a), as well as several other possibilities, e.g. thickening, slab break-off, slab drips etc (e.g. Toussaint et al., 2004a,b). Models of HP-UHP rocks exhumation can be summarized into the following groups: (1) syn-collisional

exhumation of a coherent and buoyant crustal slab, with formation of a weak zone at the entrance of the subduction channel (Chemenda et al., 1995, 1996; Toussaint et al., 2004b; Li and Gerya, 2009); (2) episodic ductile extrusion of HP-UHP rocks from the subduction channel to the surface or crustal depths (Beaumont et al., 2001; Warren et al., 2008a); (3) continuous material circulation in the rheologically weak subduction channel stabilized at the plate interface, with materials exhumed from different depths (Burov et al., 2001; Stöckhert and Gerya, 2005; Yamato et al., 2007; Warren et al., 2008a).

In this chapter, the processes and dynamics of continental subduction/collision and HP-UHP rocks exhumation are investigated by the method of large-scale numerical geodynamic modeling. First the numerical method is described, which is followed by the numerical model setup and systematic thermo-mechanical numerical experiments. The discussion section covers a broad range of topics related to the continental subduction and exhumation. Finally a concluding part is presented.

# **2. Numerical modeling method**

#### **2.1 Governing equations and numerical implementation**

The momentum, continuity and heat conservation equations for a 2D creeping flow including thermal and chemical buoyant forces are solved:

(i) Stokes equation

$$\begin{aligned} \frac{\partial \sigma'\_{\infty}}{\partial \mathbf{x}} + \frac{\partial \sigma'\_{\infty}}{\partial \mathbf{z}} &= \frac{\partial P}{\partial \mathbf{x}}\\ \frac{\partial \sigma'\_{\infty}}{\partial \mathbf{x}} + \frac{\partial \sigma'\_{\infty}}{\partial \mathbf{z}} &= \frac{\partial P}{\partial \mathbf{z}} - \mathbf{g} \,\rho(\mathbf{C}, M, P, T) \end{aligned} \tag{1}$$

where the density *ρ* depends on composition (*C*), melt fraction (*M*), pressure (*P*) and temperature (*T*); *g* is the acceleration due to gravity.

(ii) Conservation of mass is approximated by the incompressible continuity equation

$$\frac{\partial v\_{\chi}}{\partial \chi} + \frac{\partial v\_{z}}{\partial z} = 0 \tag{2}$$

(iii) Heat conservation equations

$$\rho \, \mathbf{C}\_p \left( \frac{DT}{Dt} \right) = -\frac{\partial q\_x}{\partial x} - \frac{\partial q\_z}{\partial z} + H\_r + H\_d + H\_s + H\_L \tag{3}$$

$$\begin{bmatrix} \Box & \Box & q\_x = -k(\mathbf{C}, P, T) \frac{\partial T}{\partial \mathbf{x}'} & q\_z = -k(\mathbf{C}, P, T) \frac{\partial T}{\partial z} \\\ H\_d = T a \frac{\partial P}{\partial t}, \ H\_s = \sigma\_{\text{xx}}' \dot{\varepsilon}\_{\text{xx}} + \sigma\_{\text{zz}}' \dot{\varepsilon}\_{\text{zz}} + 2 \sigma\_{\text{xx}}' \dot{\varepsilon}\_{\text{xx}} \end{bmatrix} \tag{4}$$

where *D*/*Dt* is the substantive time derivative. *x* and *z* denote the horizontal and vertical directions, respectively. The deviatoric stress tensor is defined by *σ* ′ *xx* , *σ* ′ *xz* , *σ* ′ *zz* , whilst the strain rate tensor is defined by *ε* ˙ *xx*, *ε* ˙ *xz*, *ε* ˙ *zz*. *qx* and *qz* are heat flux components. *ρ* is the density. *k*(*C*, *P*, *T*) is the thermal conductivity as a function of composition (*C*), pressure (*P*) and temperature (*T*). *C<sup>p</sup>* is the isobaric heat capacity. *Hr*, *Ha*, *Hs*, *H<sup>L</sup>* are radioactive, adiabatic, shear and latent heat production, respectively (see Table 1 for details of these parameters). To solve the above equations, the I2VIS code is used (Gerya and Yuen, 2003a). It is a two-dimensional finite difference code with marker-in-cell technique which allows for non-diffusive numerical simulation of multi-phase flow in a rectangular fully staggered

Eulerian grid. I2VIS accounts for visco-plastic deformation and several geological processes that are described below. All abbreviations and units used in this chapter are listed in Table 1.

Numerical Geodynamic Modeling of Continental Convergent Margins 3 Numerical Geodynamic Modeling of Continental Convergent Margins 275


Table 1. Abbreviations and units of the variables used in this chapter.

### **2.2 Boundary conditions**

For the 2D numerical models presented in this chapter, the velocity boundary conditions are free slip at all boundaries except the lower one, which is permeable (Burg and Gerya, 2005;

Li et al., 2010). Infinity-like external free slip conditions along the lower boundary imply free slip to be satisfied at 1000 *km* below the bottom of the model. As for the usual free slip condition, external free slip allows global conservation of mass in the computational domain and is implemented by using the following limitation for velocity components at the lower boundary: *<sup>∂</sup>vx*/*∂<sup>z</sup>* = 0, *<sup>∂</sup>vz*/*∂<sup>z</sup>* = −*vz*/∆*zexternal* , where <sup>∆</sup>*zexternal* is the vertical distance from the lower boundary to the external boundary where free slip (*∂vx*/*∂z* = 0, *v<sup>z</sup>* = 0) is satisfied. The thermal boundary conditions have a fixed value (0 ◦C) for the upper boundary and zero horizontal heat flux across the vertical boundaries. For the lower thermal boundary, an infinity-like external constant temperature condition is imposed, which allows both temperatures and vertical heat fluxes to vary along the permeable box lower boundary, implying constant temperature condition to be satisfied at the external boundary. This condition is implemented by using the limitation *<sup>∂</sup>T*/*∂<sup>z</sup>* =(*Texternal* − *<sup>T</sup>*)/∆*zexternal* where *<sup>T</sup>external* is the temperature at the external boundary and <sup>∆</sup>*zexternal* is the vertical distance from the lower boundary to the external boundary (Burg and Gerya, 2005; Li et al., 2010).

#### **2.3 Rheological model**

A viscoplastic rheology is assigned for the model in which the rheological behaviour depends on the minimum differential stress attained between the ductile and brittle fields. Ductile viscosity dependent on strain rate, pressure and temperature is defined in terms of deformation invariants as:

$$\eta\_{\text{durchile}} = (\pounds\_{II})^{\frac{1-n}{n}} F(A\_D)^{-\frac{1}{n}} \exp(\frac{E+PV}{nRT}) \tag{4}$$

where *ε* ˙*II* = 0.5 *ε* ˙*ij ε* ˙*ij* is the second invariant of the strain rate tensor. *AD*, *E*, *V* and *n* are experimentally determined flow law parameters (Table 2). *F* is a dimensionless coefficient depending on the type of experiments on which the flow law is based. For example: *F* = [2 (1−*n*)/*n* ]/[3 (1+*n*)/2*n* ] for triaxial compression and *F* = 2 (1−2*n*)/*n* for simple shear.

The ductile rheology is combined with a brittle/plastic rheology to yield an effective visco-plastic rheology. For this purpose the Mohr-Coulomb yield criterion (e.g. Ranalli, 1995) is implemented as follows:

$$
\sigma\_{yield} = \mathbb{C} + P \sin(\varphi\_{eff}) \tag{5}
$$

$$
\sin(\varphi\_{eff}) = \sin(\varphi)(1 - \lambda)
$$

$$
\overset{\prod}{\longrightarrow} \overbrace{\sigma\_{plastic}} = \frac{\sigma\_{yield}}{2\varepsilon\_{\text{II}}}
$$

where *σyield* is the yield stress. *ε* ˙ *II* is the second invariant of the strain rate tensor. *P* is the dynamic pressure. *C* is the cohesion. *ϕ* is the internal frictional angle. *λ* is the pore fluid coefficient that controls the brittle strength of fluid-containing porous or fractured media (Brace and Kohlstedt, 1980). *ϕeff* can be illustrated as the effective internal frictional angle that integrates the effects of internal frictional angle (*ϕ*) and pore fluid coefficient (*λ*). *λ* is the pore fluid coefficient that controls the brittle strength of fluid-containing porous or fractured media.

The effective viscosity of molten rocks (*M* ≥ 0.1) was calculated using the formula (Pinkerton and Stevenson, 1992; Bittner and Schmeling, 1995):

$$
\eta = \eta\_0 \exp[2.5 + (1 - M)(\frac{1 - M}{M})^{0.48}] \tag{6}
$$

where *η*<sup>0</sup> is an empirical parameter depending on rock composition, being *η*<sup>0</sup> = 10 <sup>13</sup> *Pa s* (i.e. 1 × 10 <sup>14</sup> ≤ *η* ≤ 2 × 10 <sup>15</sup> *Pa s* for 0.1 ≤ *M* ≤ 1) for molten mafic rocks and *η*<sup>0</sup> = 5 × 10 <sup>14</sup> *Pa s* (i.e. 6 × 10 <sup>15</sup> ≤ *η* ≤ 8 × 10 <sup>16</sup> *Pa s* for 0.1 ≤ *M* ≤ 1) for molten felsic rocks. Successfully tested for a broad range of suspensions with various bubble or crystal conventions, this formula takes into account, other than concentration, particle shape and size distribution.

### **2.4 Partial melting model**

The numerical code accounts for partial melting of the various lithologies by using experimentally obtained *P*-*T* dependent wet solidus and dry liquidus curves (Gerya and Yuen, 2003b). As a first approximation, volumetric melt fraction *M* is assumed to increase linearly with temperature accordingly to the following relations (Burg and Gerya, 2005):

$$M = 0, \text{ when } T \le T\_{\text{solidus}}$$

$$M = \frac{T - T\_{\text{solidus}}}{T\_{\text{liquidus}} - T\_{\text{solidus}}}, \text{ when } T\_{\text{solidus}} < T < T\_{\text{liquidus}} \tag{7}$$

$$M = 1, \text{ when } T \ge T\_{\text{liquidus}}$$

where *Tsolidus* and *Tliquidus* are the wet solidus and dry liquidus temperature of the given lithology, respectively (Table 3). Consequently, the effective density, *ρeff* , of partially molten rocks varies with the amount of melt fraction and *P*-*T* conditions according to the relations:

$$
\rho\_{eff} = \rho\_{solid} - M(\rho\_{solid} - \rho\_{molten}) \tag{8}
$$

where *ρsolid* and *ρmolten* are the densities of the solid and molten rock, respectively, which vary with pressure and temperature according to the relation:

$$
\rho\_{P,T} = \rho\_0 [1 - a(T - T\_0)][1 + \beta(P - P\_0)] \tag{9}
$$

where *ρ*<sup>0</sup> is the standard density at *P*<sup>0</sup> = 0.1 *MPa* and *T*<sup>0</sup> = 298 *K*; *α* and *β* are the thermal expansion and compressibility coefficients, respectively (Tables 1 and 3).

The effects of latent heat *H<sup>L</sup>* (e.g. Stüwe, 1995) are accounted for by an increased effective heat capacity (*CPe f f* ) and thermal expansion (*αeff* ) of the partially molten rocks (0 < *M* < 1), calculated as

*CPe f f* = *C<sup>P</sup>* + *QL*( *∂M ∂T* )*<sup>P</sup>* (10) *αeff* = *α* + *ρ Q<sup>L</sup> T* ( *∂M ∂P* )*<sup>T</sup>* (11)

where *C<sup>P</sup>* and *α* are the heat capacity and the thermal expansion of the solid crust, respectively, and *Q<sup>L</sup>* is the latent heat of melting of the crust (Table 1).

#### **2.5 Topographic model**

The spontaneous deformation of the upper surface of the lithosphere, i.e. topography, is calculated dynamically as an internal free surface by using a low viscosity (e.g. 10 <sup>18</sup> *Pa s*), initially 8-12 *km* thick layer (thickness of this layer changes dynamically during experiments) above the upper crust. The composition is either "air" (1 *kg*/*m*<sup>3</sup> , above water level) or "water" (1000 *kg*/*m*<sup>3</sup> , below water level). The interface between this weak layer and the underlying crust is treated as an internal erosion/sedimentation surface which evolves according to the

Eulerian transport equation solved in Eulerian coordinates at each time step (Gerya and Yuen, 2003b):

$$\frac{\partial z\_{\rm es}}{\partial t} = v\_z - v\_\chi \frac{\partial z\_{\rm es}}{\partial \chi} - v\_\chi + v\_\ell \tag{12}$$

where *zes* is the vertical position of the surface as a function of the horizontal distance *vx*. *vz* and *vx* are the vertical and horizontal components of the material velocity vector at the surface. *vs* and *ve* are the sedimentation and erosion rates, respectively, which correspond to the relation: *v<sup>s</sup>* = 0, *v<sup>e</sup>* = *ve*0, when *zes* < erosion level; *v<sup>s</sup>* = *vs*0, *v<sup>e</sup>* = 0, when *zes* > erosion level; where *ve*<sup>0</sup> and *vs*<sup>0</sup> are the imposed constant large scale erosion and sedimentation rates, respectively. The code allows for marker transmutation that simulates erosion (rock markers are transformed to weak layer markers) and sedimentation (weak layer markers are transformed to sediments).

# **3. Numerical model design**

Fig. 1. Initial model configuration and boundary conditions. (a) Enlargement (1700 × 670 *km*) of the numerical box (4000 × 670 *km*). Boundary conditions are indicated in yellow. (b) The zoomed domain of the subduction zone. White lines are isotherms measured in ◦C. (c) The colorgrid for different rock types, with: 1-air; 2-water; 3,4-sediment; 5-upper continental crust; 6-lower continental crust; 7-upper oceanic crust; 8-lower oceanic crust; 9-lithospheric mantle; 10-athenospheric mantle; 11-weak zone mantle; 13,14-partially molten sediment (3,4); 15,16-partially molten continental crust (5,6); 17,18-partially molten oceanic crust (7,8). The partially molten crustal rocks (13, 14, 15, 16, 17, 18) are not shown in this figure, which will appear during the evolution of the model. In our numerical models, the medium scale layering usually shares the same physical properties, with different colors used only for visualizing plate deformation. Detailed properties of different rock types are shown in Tables 2 and 3.

Large scale models (4000 × 670 *km*, Fig. 1) are designed for the study of dynamic processes from oceanic subduction to continental collision associated with HP-UHP rocks formation and exhumation. The non-uniform 699 × 134 rectangular grid is designed with a resolution varying from 2 × 2 *km* in the studied collision zone to 30 × 30 *km* far away from it. The lithological structure of the model is represented by a dense grid of 7 million active Lagrangian markers used for advecting various material properties and temperature (Gerya et al., 2008; Li et al., 2010). The subducting plate is pushed rightward by prescribing a constant


Table 2. Viscous flow laws used in the numerical experiments. (*a*) *η*<sup>0</sup> is the reference viscosity, which is calculated as *η*<sup>0</sup> =(1/*AD*) × 10 6*n* . Other parameters are illustrated in Table 1. (*b*) The molten felsic (*F* ∗ ) is used for partial molten sediment and continental crust. The molten mafic (*G* ∗ ) is used for partial molten oceanic crust. The rheological data in this table are from Ranalli (1995).


Table 3. Material properties used in the numerical experiments. (*a*) *M*<sup>1</sup> is used for sediment and continental upper crust. *M*<sup>2</sup> is used for continental lower crust. *M*<sup>3</sup> is used for oceanic crust. *M*<sup>4</sup> is used for lithospheric and athenospheric mantle. Numbers in the brackets are corresponding to material colors in Figure 1. (*b*)

*k*<sup>1</sup> =[0.64 + 807/(*T<sup>K</sup>* + 77)] exp(0.00004*PMPa*),

$$k\_2 = \left[1.18 + 474 / (T\_K + 77)\right] \exp(0.00004 P\_{MPa}),$$

*k*<sup>3</sup> =[0.73 + 1293/(*T<sup>K</sup>* + 77)] exp(0.00004*PMPa*).(*c*)

*TS*<sup>1</sup> = {889 + 17900/(*P* + 54)+ 20200/(*P* + 54) 2 *at P* < 1200 *MPa*} *or* {831 + 0.06*PatP* > 1200 *MPa*}, *TS*<sup>2</sup> = {973 − 70400/(*P* + 354)+ 778 × 10 <sup>5</sup>/(*P* + 354) 2 *at P* <

1600 *MPa*} *or* {935 + 0.0035*P* + 0.0000062*P* 2 *at P* > 1600 *MPa*}.(*d*) *TL*<sup>1</sup> = 1262 + 0.09*P*, *TL*<sup>2</sup> = 1423 + 0.105*P*.(*e*) Parameters of viscous flow laws are shown in Table 2. (*f*) This column shows the values of sin(*φeff* ), which is the effective internal frictional angle implemented for plastic rheology. The plastic cohesion is zero in all the experiments. (*g*) 1=(Turcotte and Schubert, 1982); 2=(Bittner and Schmeling, 1995); 3=(Clauser and Huenges, 1995); 4=(Ranalli, 1995); 5=(Schmidt and Poli, 1998). In this table, meanings of all the variables are shown in Table 1. Thermal expansion coefficient *α* = 3 × 10 <sup>−</sup><sup>5</sup> *K* <sup>−</sup><sup>1</sup> and Compressibility coefficient *β* = 1 × 10 <sup>−</sup><sup>5</sup> *MPa* <sup>−</sup><sup>1</sup> are used for all the rocks.

convergence velocity (*Vx*) in a small internal domain that remains fixed with respect to the Eulerian coordinate (Fig. 1).

In the numerical models, the driving mechanism of subduction is a combination of plate push (prescribed rightward convergence velocity) and increasing slab pull (temperature-induced density contrast between the subducted lithosphere and surrounding mantle). This type of boundary condition is commonly used in numerical models of subduction and collision (e.g. Toussaint et al., 2004b; Burg and Gerya, 2005; Currie et al., 2007; Yamato et al., 2007; Warren et al., 2008b) and assumes that in the globally confined three-dimensional system of plates, local external forcing coming either from different slabs or from different sections of the same laterally non-uniform slab can be significant.

Following previous numerical studies performed with similar geodynamic settings (e.g. Warren et al., 2008a; Li and Gerya, 2009) we design numerical models consisting of three major domains (from left to right, Fig. 1): (1) a pro-continental domain, (2) an oceanic domain, and (3) a retro-continental domain. The subducting pro-continent comprises a marginal unit and an interior unit. In the continental domain, the initial material field is set up by a 35 *km* thick continental crust composed of sediment (6 *km* thick), upper crust (14 *km* thick) and lower crust (15 *km* thick), overlying the lithospheric mantle (85 *km* thick) and subjacent mantle (540 *km* thick). The oceanic domain comprises an 8 *km* thick oceanic crust overlying the lithospheric mantle (82 *km* thick) and subjacent mantle (570 *km* thick). The material properties of all layers (Fig. 1) are listed in Table 3. The initial thermal structure of the lithosphere (white lines in Fig. 1) is laterally uniform with 0 ◦C at the surface and 1300 ◦C at the bottom of the lithospheric mantle (both continental and oceanic). The initial temperature gradient in the asthenospheric mantle is around 0.5 ◦C/*km*.

## **4. Model result**

#### **4.1 Reference model**

The reference model is designed with prescribed convergence velocity (*Vx*)of5 *cm*/*y*. All the other configurations and parameters are shown in Figure 1 and Table 3.

At the initial stages, the relatively strong oceanic plate subducts along the weak zone to the mantle (Fig. 2a).The continental margin subducts to >100 *km* depth, following the high-angle oceanic subduction channel (Fig. 2b). The significant characters are the detachment of subducting upper/middle crust at the entrance zone of the subduction channel with a series of thrust faults formed (Fig. 2b-d). A small amount of crustal rocks located in the lower part of the channel are detached from the plate at asthenospheric depths, indenting into the mantle wedge and forming a compositionally buoyant plume (Fig. 2c). Such sub-lithospheric plumes are discussed in detail in Currie et al. (2007) and Li and Gerya (2009). In addition, a partially molten plume forms in the deeply subducted oceanic plate and moves up vertically until it collapses at the bottom of the overriding lithospheric mantle (Fig. 2c,d). As subduction continues, another partially molten plume forms in the deeply subducted continental plate. It also moves up vertically until it collapses at the bottom of the overriding lithospheric mantle (Fig. 2e,f). The characteristics and 2D and 3D dynamics of this kind of plume are studied in detail in Gerya and Yuen (2003b) and Zhu et al. (2009).

As continental subduction continues, partially molten rocks accumulated in the subduction channel extrude upward to the crustal depths (Fig. 2d,e). Then these UHP rocks exhume buoyantly to the surface forming a dome structure (Fig. 2f). The exhumed UHP rocks are mainly located near the suture zone with a fold-thrust belt formed in the foreland extending for about 300-400 *km* (Fig. 2f). P-T paths (Fig. 2, inset) show that peak *P*-*T* conditions

Fig. 2. Enlarged domain evolution (1300 × 600 *km*) of the reference model. Colors of rock types are as in Figure 1. Time (Myr) of shortening is given in the figures. White numbered lines are isotherms in ◦C. Small colored squares indicate positions of representative markers (rock units) for which *P*-*T* paths are shown (inset). Colors of these squares are used for discrimination of marker points plotted in *P*-*T* diagrams and do not correspond to the colors of rock types.

Fig. 3. Peak metamorphic conditions of the reference model. (a) Peak pressure condition in *GPa*; (b) Peak temperature condition in ◦C.

of the exhumed rocks are 2.5-4 *GPa* and 600-800 ◦C, respectively (also see Figure 3 for the peak pressure and temperature conditions of the collision zone). This indicates the UHP metamorphic rocks are formed and exhumed from a depth >100 *km*.

# **4.2 Models with variable convergence velocity**

The reference model is further investigated with lower convergence velocity (2.5 *cm*/*y*) and higher convergence velocity (10 *cm*/*y*). All the other parameters are the same as in Tables 2

Fig. 4. Enlarged domain evolution (800 × 225 *km*) of the model with lower convergence velocity *V<sup>x</sup>* = 2.5 *cm*/*y*. Colors of rock types are as in Figure 1. Time (*Myr*) of shortening is given in the figures. White numbered lines are isotherms in ◦C. Small colored squares indicate positions of representative markers (rock units) for which *P*-*T* paths are shown (right). Colors of these squares are used for discrimination of marker points plotted in *P*-*T* diagrams and do not correspond to the colors of rock types.

In the slow convergence regime, the continental margin also subducts to >100 *km* depth along the high-angle oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 4a). The subducting upper/middle crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 4a-d). With continued continental subduction, partially molten rocks accumulated in the subduction channel extrude upward to the crustal depth (Fig. 4c,d). *P*-*T* paths (Fig. 4) show that peak *P*-*T* conditions of the exhumed rocks are 2-4 *GPa* and 600-800 ◦C.

Numerical Geodynamic Modeling of Continental Convergent Margins 11 Numerical Geodynamic Modeling of Continental Convergent Margins 283

Fig. 5. Enlarged domain evolution (1075 × 275 *km*) of the model with higher convergence velocity *V<sup>x</sup>* = 10 *cm*/*y*. Colors of rock types are as in Figure 1. Time (*Myr*) of shortening is given in the figures. White numbered lines are isotherms in ◦C. Small colored squares indicate positions of representative markers (rock units) for which *P*-*T* paths are shown (right). Colors of these squares are used for discrimination of marker points plotted in *P*-*T* diagrams and do not correspond to the colors of rock types.

In the fast convergence regime, the continental domain continues subducting along the high-angle oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 5a). Then the crustal rocks in the lower part of the channel detach from the plate at asthenospheric depths, intrude into the mantle wedge, and form a horizontal compositionally buoyant plume (Fig. 5a-c). In addition, a partially molten plume forms in the deeply subducted plate and moves up vertically until it collapses at the bottom of the overriding lithospheric mantle (Fig. 5c,d), which is similar to behavior of the reference model. After the convergence ceases (1500 *km* shorting, 15 *Myr*), the subducted continental crustal rocks in the sub-lithospheric plume extrude upward to the surface forming a dome structure (Fig. 5d). *P*-*T* paths (Fig. 5) indicate that peak *P*-*T* conditions of the exhumed rocks are 3-4 *GPa* and 600-800 ◦C, respectively. This parameter sensitivity studies indicate that the slower convergence produces very small sub-lithospheric plume (Fig. 4a), coupled subduction channel and wide collision zone (Fig. 4d). In contrast, the faster convergence results in very large sub-lithospheric plume (Fig. 5a), decoupled subduction channel and narrow collision zone (Fig. 5d). Both of the models can obtain UHP rocks exhumation. However, the convergence velocity changes the amount of crustal rocks subducted to and exhumed from UHP depth (c.f. Figs 4 and 5).

# **4.3 Models with variable thermal structure of the oceanic lithosphere**

Fig. 6. Enlarged domain evolution (900 × 225 *km*) of the model with higher temperature of the oceanic lithosphere (hot model). Colors of rock types are as in Figure 1. Time (*Myr*)of shortening is given in the figures. White numbered lines are isotherms in ◦C. Small colored squares indicate positions of representative markers (rock units) for which *P*-*T* paths are shown (right). Colors of these squares are used for discrimination of marker points plotted in *P*-*T* diagrams and do not correspond to the colors of rock types.

The lithospheric thermal structure plays an important role on subduction/collision processes (e.g. Toussaint et al., 2004a, b). Therefore we investigate the sensitivity of oceanic thermal gradient for the reference model. In the hot model, initial thermal structure of the oceanic lithosphere is linearly interpolated with 0 ◦C at the surface (≤ 10 *km* depth) and 1300 ◦Cat70 *km* depth (compared to 1300 ◦C at 100 *km* depth in the reference model). In contrast, the initial thermal structure of oceanic lithosphere in the cold model is linearly interpolated with 0 ◦Cat the surface (≤ 10 *km* depth) and 1300 ◦Cat130 *km* depth.

In the hot model, the continental margin subducts following the oceanic subduction channel to the bottom of the lithospheric mantle (Fig. 6a). Then the crustal rocks in the lower part of the channel detach from the plate at asthenospheric depths, intrude into the mantle wedge, and form a horizontal compositionally buoyant plume (Fig. 6b). The subducting upper/middle crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 6b,c). With continued continental subduction, partially molten rocks in the middle channel extrude upward to the crustal depth (Fig. 6c,d). The subduction channel is highly coupled. As a result, the partially molten rocks in the sub-lithospheric plume stay at the bottom of the overriding lithosphere (without exhumation).

Numerical Geodynamic Modeling of Continental Convergent Margins 13 Numerical Geodynamic Modeling of Continental Convergent Margins 285

Fig. 7. Enlarged domain evolution (900 × 225 *km*) of the model with lower temperature of the oceanic lithosphere (cold model). Colors of rock types are as in Figure 1. Time (*Myr*)of shortening is given in the figures. White numbered lines are isotherms in ◦C. Small colored squares indicate positions of representative markers (rock units) for which *P*-*T* paths are shown (right). Colors of these squares are used for discrimination of marker points plotted in *P*-*T* diagrams and do not correspond to the colors of rock types.

In the cold model, the continental margin subducts following the oceanic plate to the bottom of the lithospheric mantle (Fig. 7a). In this case, there is no sub-lithospheric plume formed. Consequently, the subduction channel is thicker and thicker. The subducting upper/middle crust at the entrance zone of the subduction channel detaches with thrust faults formed (Fig. 7b,c). The subducted continental crustal rocks extrude upward to the surface forming a dome structure (Fig. 7c,d). The subduction channel is highly decoupled.

The shape and characteristics of the subduction channel in the hot model (Fig. 6) is similar to that in the slow convergence model (Fig. 4). It indicates that both the higher temperature and the slower convergence can increase the rheological coupling at plate interface. As a result, coupled subduction channel is produced in these two models. In contrast, decoupled channels are formed in the colder model (Fig. 7) as well as in the faster convergence model (Fig. 5).

# **5. Discussion**

# **5.1 Flow modes in the subduction channel**

To a first approximation, viscous channel flow can be analysed using lubrication theory (e.g. England and Holland, 1979; Cloos, 1982; Cloos and Shreve, 1988a, 1988b; Mancktelow, 1995;

Raimbourg et al., 2007; Warren et al., 2008b; Beaumont et al., 2009). Under the lubrication approximations, channel flow velocity is

$$u(x,y) = -\frac{1}{2\,\eta} \frac{\partial P}{\partial x}(y\,h - y^2) + \mathcal{U}(1 - \frac{y}{h}) \tag{13}$$

where *η* is the assumed uniform viscosity in the subduction channel, *∂P*/*∂x* is the effective down-channel pressure gradient, with *x* measured in the down-dip direction, *y* is the position in the channel measured normal to the base, *h* is channel thickness and *U* is the subduction velocity of the underlying lithosphere (Fig. 8a). The overlying lithosphere is assumed to be stationary.

Fig. 8. Schematic diagram showing subduction/exhumation channel flow behavior in terms of dominating Couette (subduction) and Poiseuille (exhumation) flows (after Warren et al., 2008b; Beaumont et al., 2009). (a) General nomenclature. (b-d) Flow types identified in the models. (b) Couette flow (subduction) dominates. All flow is directed downward. (c) Buoyant materials stagnate at bottom of channel with the Poiseuille flow effect increasing. It is characterized as the transition from subduction-dominated to exhumation-dominated channel (d) Poiseuille flow (exhumation) dominates. Buoyancy-driven exhumation starts at channel bottom and propagates upward.

When nondimensional variables *u* ′ = *u*/*U*, *h* ′ = *h*/*H*, *x* ′ = *x*/*h* and *y* ′ = *y*/*h* are used, Equation 13 reduces to

$$u' = -\frac{\operatorname{E}h'^2(y'-y'^2)}{2} + (1-y')\tag{14}$$

where

$$E = \frac{H^2(\partial P/\partial x)}{\eta\_{eff}U} \tag{15}$$

$$H = (\frac{2\,\eta\_{eff}\,\mathrm{U}}{\partial P/\partial \mathrm{x}})^{\frac{1}{2}}\tag{16}$$

where *ηeff* , the effective viscosity and *∂P*/*∂x* are averaged quantities measured at the channelscale used to estimate *H*. This parameter *H* is the characteristic channel thickness for *E* ∼ 1, the balance point between downward and return flows.

Numerical models of the subduction channels can be conveniently interpreted in terms of the characteristic exhumation number *E* (Raimbourg et al., 2007; Warren et al., 2008b; Beaumont et al., 2009). The corresponding flow modes associated with burial and exhumation of UHP rocks are shown in Figure 8. The first order dynamics can be approximated by the above-mentioned lubrication theory for creeping flows, and characterized in terms of the competition between down-channel Couette flow (*U*(1 − *y*/*h*) in Eq. 13), caused by the drag of the subducting lithosphere and the opposing up-channel Poiseuille flow ((1/2 *η*)(*∂P*/*∂x*)(*yh* − *y* 2 ) in Eq. 13), driven by the buoyancy of low-density subducted crustal material. This competition can be expressed through the exhumation number *E*, which is a force ratio derived from the non-dimensional channel flow equation (Eq. 14). The actual values that determine *E* (Eqs 15, 16) will depend on the particular problem and its evolving solution. For a channel with deformable walls and no tectonic over/under-pressure, *∂P*/*∂x* =(*ρ<sup>m</sup>* − *ρc*) *g* sin(*θ*) where *θ* is channel dip (Fig. 8a).

Along with the characteristic *E*, defined at the scale of the subduction channel, space-time variations in the channel flow can be interpreted in terms of the local exhumation number *E*(*x*, *t*) and corresponding flow modes (Warren et al., 2008b; Beaumont et al., 2009). During continental subduction, *E*(*x*, *t*) evolves from <1 during subduction (c.f. Fig. 8b and Fig. 2a), to ∼1 during detachment and stagnation in the subduction channel (c.f. Fig. 8c and Fig. 2b,c), to >1 at the onset of and during exhumation (c.f. Fig. 8d and Fig. 2d-e). Buoyancy is a necessary, but not sufficient, condition for UHP exhumation. Among other controlling factors (Fig. 8), decreasing viscosity (*ηeff* ) is typically most important for driving *E* beyond the exhumation threshold. In general, *E*(*x*, *t*) should be regarded as a measure of local exhumation potential, even where the local threshold value is exceeded (*E*>1), efficient exhumation may be impeded by constrictions (small *h*) or high viscosities (large *ηeff* ) further up the channel.

#### **5.2 Coupled and decoupled subduction channel**

The numerical results show that the coupled subduction channel favors lower convergence velocity (Fig. 4) and hotter oceanic lithosphere (Fig. 6). It is characterized by continuous accretion of the weak upper continental crust resulting in the development of a thick and broad crustal wedge. In contrast, the higher convergence velocity (Fig. 5) and colder oceanic lithosphere (Fig. 7) result in decoupling of the convergent plates. Transition from coupled to decoupled regime occurs always at the early stages of continental collision indicating that insertion of rheologically weak crustal material in the subduction channel is critical for the subsequent evolution of the collision zone (Faccenda et al., 2009). The numerical models confirm that HP-UHP complexes can be formed in both coupled and decoupled channels in the wide range of convergence scenarios (Figs 2-7).

As discussed in Faccenda et al. (2009), coupled collision zones (which can be either retreating or advancing) are characterized by a thick crustal wedge and compressive stresses (i.e. Himalaya and Western Alps), while decoupled end-members (which are always retreating)

are defined by a thin crustal wedge and bi-modal distribution of stresses (i.e. compressional in the foreland and extensional in the inner part of the orogen, Northern Apennines).

**5.3 Thrust fault formation and exhumation of (U)HP units**

Fig. 9. Highly-compressional regime of continental subduction (low pull force) after Chemenda et al. (1995, 2000). (a), (b) and (e) are successive stages of continental subduction in experiments without erosion. (a)-(d) show continental subduction in experiments with erosion. 1, overriding plate; 2, upper crust; 3, lower crust, 4, eroded material (sediments).

One of the most important characteristics of the numerical models presented in this study is the formation of the thrust faults (rheological weak zones), which is followed by exhumation processes (e.g. Fig. 2). Similar detachment phenomenon is also documented in analogue models of continental subduction (Fig. 9; Chemenda et al., 1995, 1996, 2000).

The behavior of the subducted continental crust depends on two competing effects: upward buoyancy and downward subduction drag as discussed in Section 5.1. Subduction drag within the crust and mantle drives the subduction of buoyant crustal materials into larger depths (Fig 2a). At the same time the buoyancy forces and also the deviatoric stresses increase in the subduction channel. When the materials are no longer strong enough to sustain the accumulated buoyancy and deviatoric stresses, the subducted continental crust will yield with forming the rheological weak zone (thrust fault) (Fig. 2b,c) followed by the detachment and exhumation of the buoyant crustal materials (Fig. 2d-f) and release of the accumulated buoyancy and deviatoric stresses.

#### **5.4 Upper crustal structure of the HP-UHP terrane**

The upper-crustal settings of many UHP terranes share a number of structural characteristics (Fig. 10a; Beaumont et al., 2009): (1) a dome structure cored by the UHP nappe, (2) domes flanked by low-grade accretionary wedge and/or upper crustal sedimentary rocks, (3) overlying and underlying medium- to high-pressure nappes, (4) suture zone ophiolites and (5) foreland-directed thrust faults and the syn-exhumation normal faults. Our numerical models reproduce the general characteristic upper crustal structures (Fig. 10b), especially the dome structure of the HP-UHP cores, the flanked and overlaid low-grade accretionary wedge, Numerical Geodynamic Modeling of Continental Convergent Margins 17 Numerical Geodynamic Modeling of Continental Convergent Margins 289

Fig. 10. (a) General characteristics of UHP complexes and surrounding upper crust (Beaumont et al., 2009). 1, structural dome cored by UHP nappe; 2, overlying lower grade rocks; 3, suture zone ophiolites; 4, medium- to high-pressure nappes; 5, foreland-directed thrust faults; 6, syn-exhumation normal faults. Variable scale reflects size range of UHP complexes. (b) Structure of UHP units and surrounding upper crust in the reference numerical model (Fig. 2f).

the foreland-directed thrust faults and the fold-thrust belt. The ophiolites are distributed near the suture zone in the reference model (Fig. 10b).

# **5.5 Influence of tectonic overpressure on the** *P***-***T* **paths of (U)HP rocks**

The principle of lithostatic pressure is habitually used in metamorphic geology to calculate burial/exhumation depth from pressure given by geobarometry. However pressure deviation from lithostatic, i.e. tectonic overpressure/underpressure due to deviatoric stress and deformation, is an intrinsic property of flow and fracture in all materials, including rocks under geological conditions (e.g. Rutland, 1965; Brace et al., 1970; Mancktelow, 1993, 1995, 2008; Petrini and Podladchikov, 2000). Therefore, one important question is whether the principle of lithostatic pressure is applicable in subduction/collision zones where crystallization and exhumation of HP-UHP rocks take place. Some authors have argued that rocks under geological conditions are too weak to support significant overpressure (e.g. Brace et al., 1970; Ernst, 1971; Burov et al., 2001; Renner et al., 2001; Green, 2005). Yet, Stuwe ¨ and Sandiford (1994) suggested that petrologically derived P-T paths may not record depth changes only but stress changes. In addition, lithospheric-scale numerical models reveal regions where pressure may be hundreds of MPa or even several GPa higher or lower than lithostatic values (Mancktelow, 1993, 1995; Petrini and Podladchikov, 2000; Burg and Gerya, 2005). The analytical solutions show that the tectonic overpressure can be as high as 60% of the lithostatic value in the brittle regime. In contrast, the ductile overpressure is normally <10% of the lithostatic value (Li et al. 2010).

Li et al. (2010) conducted systematic numerical simulations of continental subduction/collision zones with variable brittle and ductile rheologies of the crust and mantle. In the numerical model (Figs 11, 12), the uppermost lithospheric mantle that can be considered as the wall of the subduction channel shows the largest tectonic overpressure (≥1 *GPa* and ≥50% of the lithostatic pressure). However, these overpressured zones rarely

Fig. 11. Evolution of the wedge-like subduction channel within enlarged 530 × 210 *km* domain of the original 4000 × 670 *km* model. Time (Myr) of shortening is given in the figures. White numbered lines are isotherms in ◦C. Small coloured squares (with '+' in them) show positions of representative markers (rock units) for which *P*-*T* paths are shown (right). Colours of these squares are used for discrimination of marker points plotted in the *P*-*T* diagrams and do not correspond to the colours of rock types. The solid curves show *P*-*T* paths with dynamic pressure, while the dashed curves show that with lithostatic pressure.

or never participate in the exhumation processes. Hence, they do not influence the *P*-*T* conditions of geologically distributed HP-UHP rocks in nature. The main overpressure region that may influence the *P*-*T* paths of HP-UHP rocks is located in the bottom corner of the wedge-like confined channel (Fig. 11) with the characteristic magnitude of pressure deviation on the order of ∼0.3 *GPa* and 10-20% from the lithostatic values (Fig. 12). The degree of confinement of the subduction channel is the key factor controlling this magnitude. The models show that the overpressure is small (∼10% lithostatic) and should not affect in a crucial way the metamorphic mineral equilibria of the exhumed UHP rocks. The challenge would be to identify the geological record to actually measure precisely such minor deviations. An important unresolved issue concerns the question of how the tectonic overpressure could be potentially recorded by mineral equilibria in natural rocks. Thermodynamic tools used for thermobarometry of natural rocks are mainly based on experimental data obtained under conditions of an isotropic stress state (i.e. in the absence of significant deviatoric stresses) and may not be directly applicable for recording dynamic pressure in strongly stressed rocks. Indeed, not all overpressured rocks should be strongly internally stressed. For example, weak (e.g. reacting, fluid rich) rock inclusions in strong overpressured stressed lithosphere will also have strong overpressure, but the stress state will

Fig. 12. The tectonic overpressure evolution of the wedge-like subduction model (Fig. 11). (a-d): absolute overpressure in *GPa*; (a'-d'): relative overpressure (ratio between tectonic overpressure and lithostatic pressure) in percent (%). The time (*Myr*) of shortening and tracing markers are the same as in Fig. 11.

be isotropic, i.e. similar to the conditions of laboratory experiments. This isotropic stress (dynamic pressure) will be notably different from the lithostatic pressure which will be then directly recorded by mineral equilibria of such rock inclusions. Obviously further efforts are needed to experimentally study mineral equilibria in stressed rocks.

The tectonic overpressure is also investigated for several different subduction/collision and exhumation scenarios (Fig. 13). In these numerical models, the overpressures are not significant in the mature subduction channel and/or the inner collision belt, which suggests that the overpressure that can possibly affect the HP-UHP rocks is mainly related to the wedge-like confined subduction channel (Figs 11, 12).

#### 20 Will-be-set-by-IN-TECH 292 Earth Sciences

Fig. 13. The tectonic overpressure in variable subduction-collision-exhumation scenarios.

# **6. Conclusions and future perspectives**

Continental collision was investigated with numerical models where an advancing oceanic/continental plate subducts under a fixed continent. The most important results can be summarized as follows:

(1) During the processes from continental subduction to exhumation, the flow mode in the subduction channel changes from Couette-dominated to Poiseuille-dominated flows. Numerical models of the subduction channels can be conveniently interpreted in terms of the characteristic exhumation number.

(2) The coupled subduction channel is formed in the models with slower convergence or hotter oceanic lithosphere. Whereas faster convergence or colder oceanic lithosphere result in decoupling of the converging plates.

(3) The continental margin can subduct following the oceanic plate to >100 *km* depth, and then exhume to the surface forming a HP-UHP dome. The upper crustal structures of the collision zone in the numerical models are consistent with both the analogue model and the natural UHP zones. The exhumation of UHP rocks occurs in a large variety of numerical parameters. (4) The main tectonic overpressure region that may influence the *P*-*T* paths of HP-UHP rocks is located in the bottom corner of the wedge-like confined channel with the characteristic magnitude of pressure deviation on the order of ∼ 0.3 *GPa* and 10-20% from the lithostatic values. The degree of confinement of the subduction channel is the key factor controlling this magnitude. The tectonic overpressure should not affect in a crucial way the metamorphic

mineral equilibria of the exhumed UHP rocks. The challenge would be to identify the geological record to actually measure precisely such minor deviations.

The numerical modeling method is demonstrated to be a great tool to study the geodynamic processes in the continental convergent margins. Most of the existed numerical models in the relevant topics are based on the two-dimensional regimes with the along-strike variations ignored. However, the tectonic processes remain inherently three-dimensional. So it is quite significant to conduct 3D numerical modeling to investigate the dynamics of continental collision with applicable to the natural tectonic settings (e.g. Alpine and Himalayan collision belts).

# **7. Acknowledgments**

The research leading to this work receives funding from Crystal2Plate, a FP7-funded Marie Curie Action under grant agreement number PITN-GA-2008-215353, to ZHL. N. Ribe is thanked for fruitful discussion.

## **8. References**





**Earth Sciences** Edited by Dr. Imran Ahmad Dar

ISBN 978-953-307-861-8 Hard cover, 648 pages **Publisher** InTech **Published online** 03, February, 2012 **Published in print edition** February, 2012

The studies of Earth's history and of the physical and chemical properties of the substances that make up our planet, are of great significance to our understanding both of its past and its future. The geological and other environmental processes on Earth and the composition of the planet are of vital importance in locating and harnessing its resources. This book is primarily written for research scholars, geologists, civil engineers, mining engineers, and environmentalists. Hopefully the text will be used by students, and it will continue to be of value to them throughout their subsequent professional and research careers. This does not mean to infer that the book was written solely or mainly with the student in mind. Indeed from the point of view of the researcher in Earth and Environmental Science it could be argued that this text contains more detail than he will require in his initial studies or research.

# **How to reference**

In order to correctly reference this scholarly work, feel free to copy and paste the following:

Zhonghai Li, Zhiqin Xu and Taras Gerya (2012). Numerical Geodynamic Modeling of Continental Convergent Margins, Earth Sciences, Dr. Imran Ahmad Dar (Ed.), ISBN: 978-953-307-861-8, InTech, Available from: http://www.intechopen.com/books/earth-sciences/numerical-geodynamic-modeling-of-continental-convergentmargins

# **InTech Europe**

University Campus STeP Ri Slavka Krautzeka 83/A 51000 Rijeka, Croatia Phone: +385 (51) 770 447 Fax: +385 (51) 686 166 www.intechopen.com

# **InTech China**

Unit 405, Office Block, Hotel Equatorial Shanghai No.65, Yan An Road (West), Shanghai, 200040, China Phone: +86-21-62489820 Fax: +86-21-62489821

© 2012 The Author(s). Licensee IntechOpen. This is an open access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.